female strobili model examination

library(rstan)

## Loading required package: ggplot2

## Loading required package: StanHeaders

## rstan (Version 2.18.1, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

library(bayesplot)

## This is bayesplot version 1.5.0

## - Plotting theme set to bayesplot::theme_default()

## - Online documentation at mc-stan.org/bayesplot

library(ggplot2)
library(rethinking)

## Loading required package: parallel

## rethinking (Version 1.59)

fitfem.stan <- readRDS("2019-04-17_fem_intercept_only.rds")
femsum <- precis(fitfem.stan, depth = 2)
param_names <- fitfem.stan@model_pars
post <- as.array(fitfem.stan)
postsamples <- extract.samples(fitfem.stan)

lp <- log_posterior(fitfem.stan) #logposterior
np <- nuts_params(fitfem.stan) #nuts parameters

#divergent transitions

color_scheme_set("darkgray")

mcmc_parcoord(post, np=np, regex_pars = param_names[1:10])

mcmc_pairs(post, np=np, regex_pars = c("beta", "kappa"))

mcmc_pairs(post, np=np, regex_pars = c("prov"))

mcmc_pairs(post, np=np, regex_pars = c("site"))

mcmc_pairs(post, np=np, pars = c("a_prov[1]", "a_prov[2]", "a_site[1]", "a_site[2]" ))

color_scheme_set("mix-brightblue-gray")
mcmc_trace(post, regex_pars="kappa")

color_scheme_set("red")
mcmc_nuts_divergence(np, lp)

#energy
mcmc_nuts_energy(np)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#rhat
color_scheme_set("brightblue")
rhats <- rhat(fitfem.stan)
mcmc_rhat(rhats)

mcmc_rhat_hist(rhats)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#effective sample size
#ratios <- neff_ratio(postsamples)
#mcmc_neff(ratios, regex_pars="sigma")

#autocorrelation
mcmc_acf(post, regex_pars = "kappa", lags=10)

# lots of autocorrelation :( :())
mcmc_acf(post, regex_pars = "sigma")

mcmc_intervals(post, regex_pars = c("sigma")) + ggtitle("sigmas")

mcmc_areas(post, regex_pars="sigma") + ggtitle("sigmas")

mcmc_intervals(post, regex_pars = c("a_prov", "a_site", "a_year")) + ggtitle("intercepts")

mcmc_intervals(post, regex_pars = c("a_clone")) + ggtitle("clone intercepts")

mcmc_areas(post, regex_pars = "kappa") + ggtitle("cutpoints")

mcmc_areas(post, pars="beta") + ggtitle("transition speed")